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I. Introduction 


The following report summarizes the research carried out during the 
recently completed phase of NASA grant NCCI-99. The work described here was 
carried out by the principal investigator, James Schwing, and a graduate 
research assistant, Jan Spangler, in conjunction with the SMART system design 
team (Solid Modeling Aerospace Research Tool) of the Vehicle Analysis Branch at 
NASA Langley. 

The major effort of the past six months has been the development of 
software used with the derivation of smooth 3-D surfaces from a sequence of 
cross-sections. Additional work has considered on problems arising in the 
creation of surfaces by extrusion and the presentation of calculated physical 
properties. 



II. Minimizing Error in Surface Representation. 


A. The Problem. 

The basis for geometric representation used in the SMART system is the 
bicubic parameterization known as the Bezier patch. Refer to Foley and 
Van Dam [1] for example for a discussion of basic patch definition and manipu- 
lation. The SMART system attempts to provide interfaces for the design 
engineer that corresponds to natural engineering design and development tools. 
Thus one of the geometric input techniques provided allows the designer to 
input a sequence of cross-sections of the object under consideration. The 
problem then becomes one of converting this sequence of cross-sections to a 
collection of Bezier patches that reproduces the given data as accurately as 
possible. 

The key to this conversion lies in the calculation of Bezier curves, edges 
of the Bezier patches, which approximate the cross-sections in a way that 
minimizes error. A straight forward solution to this would seem to be the 
calculation of these Bezier curves via a least squares technique. However, the 
process is complicated by the fact that the calculated Bezier curves are then 
used as the edges of the Bezier patches that are expected to join together in a 
smooth, differentiable fashion. This imposes additional interaction conditions 
on the calculation of the least square Bezier curves. 

Description of these additional conditions can best be seen by considering 
what happens when two patches are to be joined in a differentiable fashion. 

The reader is again referred to Foley and Van Dam [ 1 ] for more detail . To 
summarize let Patch 1 be represented by the 16 control points P^ ... P 44 and 
Patch 2 be represented by Qn ... Q 44 . The conditions can now be described in 
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terms of these control points. Suppose that the common edge is given by the 
control points P14, P24> P34, P44 in Patch 1 and by the control points Qn, 
Q21, Q31, Q41 in Patch 2. The first requirement is that P^4 - Q^i for i » 1, 
2, 3, and 4. In addition, the following relation must hold for the interior 
control points: Qj.2 ■ Qii + k (Q^i - Pjj) with k constant and i » 1, 2, 3, 

and 4. 

As an example let the two curves shown below in figure 1 represent a 
portion of two consecutive cross-sections where each portion is represented by 




the joining of two Bezier curve segments. The condition cited above requires 
that 

P5 ■ P4 + k (P4 - P3) and 

Q 5 - Q4 + k (Q4 - Q3) for some k. 

Since each cross-section is adjacent to at least one other cross-section the 
value of k must be determined for all cross-sections simultaneously. 


B. The Theory. 

Assume that the Bezier curve segments which approximate a given cross- 
section are ordered and that the Bezier control points for the previous 
segments have already been determined. That is, in each cross-section there is 
a previously determined Bezier curve segment, such as P^, ... , P4, to which a 
new Bezier curve segment, such as P4, ... , P7, will be joined. The special 
case of curve segments where there is no previous information will be treated 
separately. 

Figure 2 below represents the segment to be calculated taken from say the 
jth cross-section. The given data points represent points of the jth cross- 
section. Interpolation conditions for the Bezier curve require that the first 
and fourth control points of this portion correspond to the first and last data 
points. Thus the values for Pqj and P^j are easily determined. By assumption 
stated above the control points of the previous segment are known, specifically 
cjj. The requirements noted in the previous section state that the control 
* 
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point CQj must satisfy with constant k: 

Cqj ■ P 0 j + k (P 0j - cxj) for all cross-sections j. 


Figure 2 



The least squares technique is now applied to find those values of k and 
Cij that minimize the error made by approximating the data points by the Bezier 
curves* To solve the least squares problem this portion of the curve is 
parameterized by chord length along the cross-section data points. Let those 
distances be represented by s^j. Further let lj, Lqj and L]j represent the 
distances and dj and Dj represent the direction vectors necessary to define the 
% 
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control points as follows: 


c lj - p 0j - Ij «*j 
°oj - p 0j + kij d j 
c lj “ p lj + L lj Dj 

and where it is known that Lqj « k lj. 

The least square solution then leads to both of the following types of 
equations. 


Eqn. 1: 
Eqn. 2: 


where 


aj k + bj Lij - Cj 

N N 

k + 6j Li ^ " f j 

n-i-1 


a-,- ■ 


i 

3 lj [ £_j (1 - 3 ij) 2 s ?j] [^xOj ^xl j + dyOj dylj] 

3 [^ZU - s ij) 2 s*j] t<tllj + dylj] 
i*l 

n-j-l 

S {(1 - s ij) s ?j ( x ij d xlj + yij d yljl) 
i»l 

n-,— 1 

“ (g (1 - s ij) 3 S 2j (1 + 2 SijOJtxQj d xl j + y 0 j d ylj ] 
nj-1 

C 1 ~ s ij) s xj (3-2 s ij) } [ x nj j d xlj + y n -jj d yl j3 

nj-1 

3 lj 123 (1 - Sij) 4 s ijHd3£ 0 j + d yOj J 
i=*l 
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nj-1 

lj { 2Z (1 ~ s ij)^ s ij t x ij d x Oj + yij ^yOjl 
i-1 

nj-1 

—[ (1 ~ s ij)^ s ij (1 + 2 Sj_j)][xQj d x Qj + yoj dyOj] 

n-j-1 

“ s ij)^ s ?j (3-2 Sj^j)] ( x njj d xOj + 7njj ^yOj] ) 


Eqn.s 1 and 2 can be solved for all j by the following: 

N N 

k - S[ f j - ( a j' c j)/bj] / 5Z?Uj - a j/ b j] 
i»l i*l 

1*1 j “ [cj — (sj* k) ] / bj . 

For the special case that considers the first segment of each cross- 
section, the interaction condition weakens. This is a consequence of the fact 
that these patches will have no continuity condition on their leading edge. 

Thus the ratio constraint previously described no longer applies and the values 
for L may be calculated at both ends of the curve. 


Fortunately, the solution of the resulting least squares problem for this 
case leads to virtually the same coefficients, aj, bj, cj, dj, ej, fj as listed 
above. One need only set lj to 1. The solution to the equations in this 
special case is given by: 


LOj - £ b j* f j " a j‘ c jl / f b j- d j " a i3 

L lj = [cj-dj - aj-fj] / [bj-dj - aj] 
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C. The Results. 


As implemented these routines will take as input a collection of cross- 
sections and produce a best least squares approximant satisfying the conditions 
necessary for the building of a smooth Bezier surface. The appendix includes 
the code used for implementing these ideas. 


III. Solids Via Extrusion. 


Here the basic idea is that a given solid may be defined by the act of 
dragging a fixed, user-defined cross-section along another user-defined path. 

No major restrictions are placed upon either the cross-section, which may be 
any planar curve, or the path, which may be any 3-D curve. Both are repre- 
sented internally by the standard SMART format of connected Bezier curve 
segments. As in the previous section these Bezier curves are to be used to 
generate the corresponding surface patches of the solid being defined. This 
research addressed the problem of providing the same continuity and shape in 
the resulting solid as that of the underlying extrusion path. 

There are two aspects to the problem mentioned above. First, when 
segments of the underlying extrusion curve join in a continuous fashion, it is 
necessary to join the resulting surfaces in a continuous fashion. Secondly, 
the shape of the "tube” generated by the extrusion should accurately reflect 
the user-defined cross-section. The solution to the first half of this problem 
i$ discussed in the latter part of this chapter. 
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With respect to the second half of the problem, the solution is not 
immediately obvious. Unfortunately, the process is not "well defined." That 
is, it is not possible to completely determine all of the final parameters from 
the two curves described above as the user input. In effect, there is one 
remaining parameter left to be freely picked by the software. 

To this point all existing automatic techniques that have been employed to 
determine this final parameter lead to an undesirable twist in the extruded 
surface in some cases. That is, the resulting surface appears to twist so that 
the inside is totally constricted. Figure 3 below illustrates this condition. 
Currently we are still working on a solution to this problem. 



* 


Figure 3 
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The key to the first half of the problem again rests in the proper 
determination of the constant ratio k mentioned in the previous section. The 
solution proposed below not only will keep the continuity of the joining 
surfaces the same as that of the underlying extrusion curve, it also attempts 
to reproduce as faithfully as possible the shape of the underlying extrusion 
curve. In order to do that the following aesthetic was adopted: 

Aesthetic: If the underlying extrusion curve is either linear or 

circular reproduce the result exactly. 

Note first that Bezier patch control points are known at the join, since 
they are given precisely by the user-defined cross-section placed at the join 
and oriented so that it is normal to the extrusion curve. Thus the solution to 
this problem reduces to identifying the appropriate values of the patch control 
points immediately preceding and following this join. 

Let the extrusion curve preceding the join have the control points Cj, C2, 
C3, C4. Similarly let C4, C5, Cg, C7 represent the control points of the curve 
segment following the join. Define and M2 as the "midpoints" of these 
respective curve segments. 

Mi - (C x + 3 C 2 + 3 C 3 + C4) / 8 
M 2 . - (C4 + 3 C 5 + 3 C 6 + C 7 ) / 8 

Finally define W to be the center of the circle containing the three points Mj, 
M2, C4. Note that if these point are collinear, then W cannot be defined but 
that we have the simple case of reproducing the new control points in a linear 
fashion. If on the other hand the points are not collinear, then the center W 
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is used to compute an appropriate radius of curvature for each of the given 
control points, P, in the user-defined cross-section. This radius is then used 
to produce the patch control points preceding and following the join the 
correspond to P. Specifically, we determine the previously mentioned constant 
Ic from the underlying curve as follows: 

Eqn. 3 k - 1 1 C \ \ / \ \ C^Ci 1 1 . 

By using this k in the determination of all interior preceding and following 
Bezier points, the appropriate continuity class is assured. 

Let Pp and Pf be the interior control points that precede and follow P in 
the definition of the Bezier surfaces meeting at the join. Notice that the 
following relations must hold: 

Eqn. 4 Pp - P + rj (C 3 - C 4 ) 

Eqn. 5 Pf ■ P + r 2 (C 5 - C 4 ). 

As mentioned above, it is desired for continuity sake that these points share 
the common k value, that is: 

k- II ft; || / II FT? II- 

The definitions of Pp and Pf above show that 

k - ri || CTO; II / r 2 1 1 C^Ci 1 1 

- (ri / r 2 ) (|| CTC 5 || / || cTci ||) 

- (ri / r 2 ) k. 

This implies that r^ - r 2 . At this point it is possible to combine this 
requirement and Eqn.s 3, 4, and 5 with the previously mention aesthetic to 
derive values for Pp and Pf. Basically the center of curvature W is used with 
these facts via: 

• Linear Case: r i ** r 2 ■ 1 

Non-linear Case: r l a r 2 * 1 1 ^ W 1 1 / | | C 4 W | | . 
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Equations for the calculation of W are straight forward in both two and 
three dimensions and can be found in any calculus text. 


IV. Display of Physical Properties. 

Calculation of physical properties is an important step in the analysis of 
aerospace vehicles. During the process of conceptual design such calculations 
must be relatively efficient without sacrificing significant accuracy so that a 
multiplicity of ideas can be tried rapidly with a reasonable confidence in the 
results. Previous research under this grant produced mathematical software 
provided accuracy beyond the required tolerances and which proved to be up to 
four times faster. Work over the last six months on this topic involved 
developing an appropriate user interface for the display of this information. 

It has turned out that it is natural to represent the geometry for 
aerospace systems in tree data structures. These trees capture the hier- 
archical relation between system components and their subassemblies. Since the 
physical properties are calculated for each of the basic subassemblies and then 
propagated through the hierarchy to the more complex components. This infor- 
mation combined with the fact that the designer is already interacting with 
this hierarchy through the mouse dictated the following design of the inter- 
face. 

Once the designer selects the calculation of physical properties, a dual 
viewport display is presented. One of the viewports contains a representation 
of the current hierarchy. From this viewport the user may use the mouse to 
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move around the data structure and select a particular component or subassembly 
for which the properties should be presented. The selected information is then 
immediately displayed in the other viewport. In a sense this allows the 
designer to browse the model each time design changes have been made. 



REFERENCE 


Foley, J. and A. Van Dam, Principles of Interactive Computer Graphics . 
Addison Wesley, 1982. 
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ROUTINES USED IN THE GENERATION OF MINIMUM ERROR BEZIER 
CURVES CORRESPONDING TO A SET OF CROSS SECTIONS. 


34 

gen_bz 

111 

seg_info 

151 

get_dvect 

203 

bz_minerrl 

267 

ln_vect 

304 

setscl 

337 

extrap 

373 

tancalc 

4*10 

sep_seg 

452 

solvel 

495 

solve2 

542 

get_cpts 


INCLUDE FILE DEPENDENCIES 

<sgimath.h> 

. ./include. dir/act_data.h 


DEFINED FUNCTONS AND VARIABLES 

sqr(X) (X) * (X) 

MAXCPCSP1 11 /* represents MAXCPCS + 1 */ 

MAXCPCSP2 12 /* represents MAXCPCS + 2 */ 


ROUTINE: gen_bz 

PURPOSE: driving routine for the generation of minimum error 

Bezier control points for Bezier curves approximating 
a set of cross sections 

CALLING PROCEDURE/ 

DECLARATIONS : 

gen_bz (xp,yp, zp, ifl ,ncs,npcs, opts, err) 

float xp [ ] [MAXCS ] , yp [ ] [MAXCS ] , zp [ ] , *cpts [ ] [MAXCS ] [4] [3] 
int ifl [] [MAXCS] ,ity,ncs,npcs[] ,*err; 

INPUT VARIABLES: 

xp,yp - [MAXCPCS] [MAXCS] , arrays containing the cross 
section curves 

zp - [MAXCS], array containing the z position of the 
given cross section 

ifl - [MAXPPCS] [MAXCS] , array containing an indicator 
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for curve segment break points and continuity 
conditions 

ncs - actual number of cross sections used 
npcs - [MAXCS] , array listing the actual number of points 
used in each cross section 

OUTPUT VARIABLES: 

cpts - [MAXCPCS] [MAXCS] [4] [2] , array containing the x,y 
coordinates of the 4 Bezier control points for each 
err - if non-zero indicates an error exit condition 
due to inconsistent data 
curve and each cross section 

EXTERNAL VARIABLES : none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

n - [MAXCPCSP2] [MAXCS] , array containing the actual numbei 
of points in each curve segment 
nseg - the actual number of curve segments for any cross 
section 

cs, segno - counters 

brkpt - [MAXCPCSP1] [MAXCS] , array containing the index 
to xp,yp for the break points between the Bezier 
curve segments 

dvect_in , dvect_out - [MAXCPCSP1] [MAXCS] [2] , array containi 
the x,y components of the incomming/ outgoing 
tangent directions at each break point 
x,y - [MAXPPCS] [MAXCS] , array containing the coordinates c 
a selected curve segment for each cross section 
dO,dl - [MAXCS] [2], array containing the x,y components 
for the start/end derivative direction on a given 
curve segment over all cross sections 
lambda - [MAXCS] [[2 ] , array containing the length of each 
of hte control points from the start/end points 
of a given curve segment over all cross sections 
lambda_old - [MAXCS], array containing the length of the 
end control point of the previous section stored 
over all cross sections 


ROUTINES INVOKED: 

seg_info 

get_dvect 

sep_seg 

bz_minerrl 

* get_cpts 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 
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REVISIONS : 

DATES : 


ROUTINE: seg_info 

PURPOSE: routine used to obtain information about the curve 

segments that will be used to approximate the cross 
sections 

CALLING PROCEDURE/ 

DECLARATIONS * 

seg_info(ifl,ncs,npcs,n, &nseg,brkpt, &err) 

int ifl[] [MAXCS] ,ncs,npcs[ ] ,n[] [MAXCS] , *nseg,brkpt[ ] [MAXC 
*err; 

INPUT VARIABLES: 

ifl,ncs,npcs - as described in "gen_bz H above 
OUTPUT VARIABLES: 

n,nseg,brkpt,err - as described in "gen_bz H above 

EXTERNAL VARIABLES : none 

GLOBAL VARIABLES: none 

INTERNAL VARIABLES: 

cs,seg,i,j - counters 

segcnt - track the number of curve segments from cross 
section to cross section 

ROUTINES INVOKED: none 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS: 

DATES : 


ROUTINE : get_dvect 

PURPOSE: routine used to produce incoming and outgoing tangent 

directions for each curve segment break point 


CALLING PROCEDURE/ 
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DECLARATIONS Z 

get_dvect ( xp , yp , ncs , nseg , n , if 1 , brkpt , dvect_in , dv ect_out ) 

float xp [ ] [MAXCS] , yp [ ] [MAXCS ] , 

dvect in [ ] [MAXCS ] [ 2 ] , dvect_out [ ] [MAXCS ] [ 2 ] ; 
int ncs, nseg, n[] [MAXCS] ,ifl[] [MAXCS] , brkpt [] [MAXCS] ; 

INPUT VARIABLES: 

xp, yp, ncs, nseg, n, if 1, brkpt - as described in M gen_bz' 
OUTPUT VARIABLES: 

dvect_in , dvect_out - as described in "gen_bz" 

EXTERNAL VARIABLES : none 

GLOBAL VARAIBLES : none 

INTERNAL VARIABLES: 

a,b - [4], array containing tangent approximation info 
ct,st - x,y components of the tangent vector 
v - [ 2 ] , temporay vector 
len - vector length 

xel,xe2,yel,ye2 - extrapolated points 
cs , seg , i - counters 

nl — the number of points used to define the incoming 
curve segment at a given break point 
n2 - the n umb er of points used to define the outgoing 
curve segment at a given break point 
ibase - index to xp,yp for the current break point 

ROUTINES INVOKED: 

tancalc 

ln_vect 

extrap 

setscl 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS : 

DATES : 


* ROUTINE : bz_minerrl 

PURPOSE: routine used to calculate the length of the Bezier 

control points along their tangent vectors so that 
a minimum least square error for a given curve segment 
over all cross sections between the approximating 
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Bezier curve and the cross section data is obtained 

CALLING PROCEDURE/ 

DECLARATIONS : 

bz_minerr 1 ( x , y , do , dl , ncs , n , segno , lambda_old , lambda ) ; 

float x[] [MAXCS] ,y[ ] [MAXCS] ,d0[][2],dl[][2], lamda_old[ ] , *lambda[ ] [2 ] 
int ncs , n [ ] [MAXCS ] , segno ; 

INPUT VARIABLES: 

x, y, dO ,dl, ncs ,n, segno, lambda_old - as described in "gen_bz" 

OUTPUT VARIABLES: 

lambda - as described in "gen_bz M 

EXTERNAL VARIABLES : none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

dist - chord length for this segment, 
s - array for chord length position of data defining 
this segment. 

the second index represents powers of s. 
scomp - (1 - s) 

the second index represents powers. 
sx2m3 -(3-2*s) 
sx2pl - (2 * s + 1) 

11, 12, 13 - accumulators for the LHS of the matrix 
equation. 

rl, r2, r3, r4, r5, r6 - accumulators for the RHS 
of te matrix equation. 

a, b, c, d, f - matrix entries used to find the 
least squares solution. 

r - common ratio reqired for surface continuity, 

found in the solution of the matrix system, 
totl, tot2 - used in solution of the matrix system. 
n_cs - number of points for this segment and this . 

cross section, 
nml - (n_cs - 1) 
nm2 - (n_cs - 2) 

ROUTINES INVOKED: 

sqr 
sqrt 
solve 1 

, solve2 

AUTHOR: James Schwing, Old Dominion University 


DATE 


12-24-86 
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REVISIONS: 

DATES : 


ROUTINE: ln_vect 

PURPOSE: routine used to find the length of a 2-D vector 

CALLING PROCEDURE/ 

DECLARATIONS : 

ln_vect(v) 
float v[2] ; 

INPUT VARIABLES: 

v - [2], vector for which the length will be calculated 
OUTPUT VARIABLES: 

ln_vect - function value, the length 

EXTERNAL VARIALBLES : none 

GLOBAL VARIABLES : none 

INTERNAL VARAIBLES : 

len - temporary storage for length 

ROUTINES INVOKED: 

sgr 

sqrt 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS : 

DATES: 


ROUTINE: setscl 

PURPOSE: routine which modifies start/ end tangent derivatives 

'calling PROCEDURE/ 

DECLARATIONS: 

setscl (&ct,&st) 
float *ct,*st; 
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INPUT VARIABLES: none 

OUTPUT VARIABLES: 

ct,st - x,y components of the tangent vector 
EXTERNAL VARIABLES: none 

GLOBAL VARIABLES i none 

INTERNAL VARIABLES: none 

ROUTINES INVOKED: none 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS: 

DATES : 


ROUTINE : extrap 

PURPOSE: routine used to extrapolate points of a curve segment 

beyond (or prior to) its end points 

CALLING PROCEDURE/ 

DECLARATIONS: 

float extrap ( a, b r c) 
float a,b,c; 

INPUT VARIABLES: 

a,b,c - input coordinates from the curve segment; x or y 
OUTPUT VARIABLES: 

extrap - function value, projected coordinate value 

EXTERNAL VARIABLES: none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

d - temporary storage 

ROUTINES INVOKED: none 


AUTHOR 


James Schwing, Old Dominion University 




Jan 16 12 s 48 1987 xsect_to_bez .d Page 8 

DATE: 12-24-86 

REVISIONS : 

DATES : 


ROUTINE: tancalc 

PURPOSE: routine used to calculate the tangent at the break point: 

of curve segments in the cross sections 


CALLING PROCEDURE/ 
DECLARATIONS : 


tancalc (a,b, &dcos,&dsin) 
float a[4] ,b[4] ,*dcos,*dsin; 


INPUT VARIABLES: 

a,b - [4 ] , differences of x,y curve values about the 
break point 


OUTPUT VARIABLES: 

dcos,dsin - x,y components of the tangent vector 

EXTERNAL VARIABLES : none 

GLOBAL VARIABLES : none 

ROUTINES INVOKED: 

sqr 

sqrt 

fabs 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS: 

DATES S 


ROUTINE : sep_seg 

PURPOSE: routine used to separate all information concerning 

an indicated curve segment in the appropriate arrays 
for minimum error processing 


CALLING PROCEDURE/ 
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DECLARATIONS: 

sep_seg (xp,yp, dvect_in , dvect_out , brkpt , segno , ncs ,x,y,dO,dl) 

float xp[] [MAXCS] ,yp[] [MAXCS] ,dvect_in[] [MAXCS] [2] ,dvect out[ ] [MAXCS] [2] 
*X[ ] [MAXCS] , *y [ ] [MAXCS] ,d0 [ ] [2 ] , dl[ ] [2] ? 

Int brkpt [ ] [MAXCS ] , segno , ncs ; 

INPUT VARIABLES: 

xp , yp , dvect_in , dvect_out , brkpt , segno , ncs - as described 
in M gen_bz" 

OUTPUT VARIABLES: 

x,y,dO,dl - as described in M gen_bz M 

EXTERNAL VARIABLES : none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

cs,i - counters 

i start , iend - markers for the start/ end subscripts of 
the current curve segment in xp,yp 

ROUTINES INVOKED: none 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS : 

DATES : 


ROUTINE : solvel 

PURPOSE: routine used to solve the least squares equations without 

constraint for the values representing the length of the 
Bezier control points from their respective end points 

CALLING PROCEDURE/ 

DECLARATIONS: 

solvei(a,b,c,d, f ,ncs,n,x,y, lambda) 

float a[ ], b[],c[],d[],f[],x[] [MAXCS] ,y[] [MAXCS] , *lambda[ ] [2 ] ; 

* int ncs,n[] [MAXCS] ; 

INPUT VARIABLES: 

a, b, c,d,f - as described in "bz_minerrl" 
ncs,n,x,y - as described in "gen_bz" 


Jan 16 12:48 1987 xsect_to_bez .d Page 10 


OUTPUT VARIABLES: 

lambda - as described in "gen_bz" 


EXTERNAL VARIABLES : none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

cs - counter 

det - partial determinant of the least squares matrix 
v - temporary vector 

ROUTINES INVOKED: 

sqr 

ln_vect 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS: 


DATES: 


ROUTINE: SOlve2 

PRUPOSE: routine used to solve the least squares equations 

constrained so that "lambda [0] / lambda__old" is 
constant over all cross sections; first solving for 
"r" the appropriate value of that ratio then "lambda" 
the values representing the length of the Bezier control 
points from their respective end points 

CALLING PROCEDURE/ 

DECLARATIONS : 

solve2 ( a , b , c , d , f , ncs , segno , n , lambda_old , x , y , lambda ) 

float a[] ,b[] /C[] ,d[] ,f [] ,lambda_old[] ,x[] [MAXCS] ,y[] [MAXCS] , 
♦lambda [] [2] ; 

int ncs , segno , n [ ] [MAXCS ] ; 

INPUT VARIABLES: 

a,b,c,d, f - as described in "bz_minerrl" 

ncs, segno ,n, lambda_old,x,y - as described in "gen_bz" 

OUTPUT VARIABLES: 

lambda - as described in "gen_bz" 


EXTERNAL VARIABLES: 


none 


Jan 16 12:48, 1987 xsect_to_bez . d Page 11 


GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

totl,tot2 - accumulators 
r - ratio described above 
cs,segpl - counters 
v - temproary vector 

ROUTINES INVOKED: 

ln_vect 

AUTHOR: James Schwing, Old Dominion University 

DATE: 12-24-86 

REVISIONS: 

DATES: 


ROUTINE : get_cpts 

PURPOSE: routine which converts the "lambda" lengths of Bezier 

control points to coordinates for a given curve 
segment over all cross sections 

CALLING PROCEDURE/ 

DECLARATIONS : 

get_cpts ( xp , yp , zp , do , dl , brkpt , segno , ncs , lambda , cpt s ) 

float xp[] [MAXCS] ,yp[] [MAXCS] ,zp[] ,d0[] [2] ,dl[] [2] , lambda [] [2] 
*CptS[] [MAXCS] [4] [3] ; 
int brkpt [ ] [MAXCS ] , segno , ncs ; 


INPUT VARIABLES: 

• xp, yp , do, dl, brkpt, segno, ncs, lambda - as described in 

"gen_bz" 

OUTPUT VARIABLES : 

cpts - as described in "gen_bz" 

EXTERNAL VARIABLES: none 

GLOBAL VARIABLES : none 

INTERNAL VARIABLES: 

* cs , i - counters 

istart,iend - indecies for the start/end points of a 
given curve segment 


ROUTINES INVOKED 


none 
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AUTHOR: 

DATE: 

REVISIONS 


48 1987 xsect_to_bez.d Page 12 

James Schwing, Old Dominion University 
12-24-86 
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/* 

ROUTINES USED IN THE GENERATION OF MINIMUM ERROR BEZIER 
CURVES CORRESPONDING TO A SET OF CROSS SECTIONS. 


46 

gen bz 


95 

seg info 


145 

get dvect 


319 

bz minerrl 


424 

In vect 


440 

setscl 


460 

extrap 


476 

tancalc 


511 

sep_seg 


551 

solvel 


584 

solve2 


639 

get_cpts 

*/ 

*/ 


# include <sgimath.h> 

# include <gl.h> 

# include " . . / include . dir/act_data . h M 

#define sqr(X) (X) * (X) 

#def ine MAXCPCSP1 11 
# define MAXCPCSP2 12 

/* */ 

seg_info() ; • 
get_dvect() ; 
sep_seg ( ) ; 
bz jminerrl ( ) ; 
get_cpts() ; 
solvel ( ) ; 
solve2 ( ) ; 
float extrap () ; 
tancalc ( ) ; 
setscl ( ) ; 
float ln_vect(); 

/* */ 

gen_bz ( xp , yp , zp , if 1 , ncs , npcs , cpts , err) 

float Xp [ ] [MAXCS ] , yp [ ] [MAXCS ] , zp [ ] , cpts [ ] [MAXCS ] [4] [3] ; 
int if 1 [ ] [MAXCS ] , ncs , npcs [ ] , *err ; 

{ int n[MAXCPCSPl] [MAXCS] ,nseg,brkpt[MAXCPCSPl] [MAXCS] ,cs, segno; 
float dvect_in[MAXCPCSPl] [MAXCS] [2] ,dvect_out[MAXCPCSPl] [MAXCS] [2] , 
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X[MAXPPCS] [MAXCS] ,y[MAXPPCS] [MAXCS] ,dO [MAXCS] [2] ,dl [MAXCS] [2] , 
lambda_old [MAXCS ] , lambda [MAXCS ] [2] ; 

M 

/* separate basic curve information */ 
seg_info(ifl,ncs,npcs,n, &nseg, brkpt, err) ; 

if (err l* 0) { 

printf (" Inconsistent data; Bezier calculation terminated\n" ) ; 
return (*err) ? 

} 


/* get tangent vectors for ALL curve */ 
/* segment end points */ 

get_dvect ( xp , yp , ncs , ns eg , n, i f 1 , brkpt , dvect_in , dvect_out ) ; 

for (cs - 0 ; cs < ncs ; cs++) 
lambda_old[cs] - 1? 


/* loop over each curve segment */ 
for (segno » 0 ; segno < nseg ; segno++) { 

/* separate the curve segment info */ 
sep_seg ( xp , yp , dvect_in , dvect_out , brkpt , segno , ncs ,x,y,dO,dl) ? 

/* minimize the approximation error */ 
bz_minerr 1 ( x , y , do , dl , ncs , n , segno , lambda_old , lambda ) ; 

for (cs » 0 ; cs < ncs ; cs++) 

lambda_old[cs] * lambda [cs] [1] ; 

/* save the new control points */ 
get_cpts (xp , yp , zp , do , dl , brkpt , segno , ncs , lambda , cpts ) ; 

} /* end loop over curve segments */ 


} 


/* V 

seg_inf o ( if 1 , ncs , npcs , n , nseg , brkpt , err) 

int ifl[ ] [MAXCS] , ncs, npcs [] ,n[] [MAXCS] , *nseg, brkpt [ ] [MAXCS] , *err; 

X int cs, seg,segcnt, i, j ; 

*err * 0; 

/* loop over all cross sections */ 

for (cs » 0 ; cs < ncs ; cs++) { 
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seg - 0; 
j - 0; 


/* loop over all point of a given */ 

/* cross section */ 

for (i ■ 0 ; i < npcs[cs] ; i++) { 

/* a break point between curve */ 

/* segments is ID'd; save that info */ 

if (ifl[i] [cs] !- 0) { 
brkpt [ seg ] [cs] * i? 
n[seg] [cs] - j ; 
seg++ ; 
j * 2; 

} 

else /* not a segment end point */ 

j++; 

) 

n[seg] [cs] - 0; 
segcnt * seg - 1; 

if (cs — 0) /* insure that each cross section */ 

*nseg ■ segcnt; /* has the same number of segments */ 

else if (*nseg l** segcnt) { 

*err = -1; 

printf("* ERROR: cross section #%d has a different number of curve\n" 

cs) ; 

printf( M segments than cross section #l\n") ; 

) /* end loop over points */ 

} /* end loop over cross sections */ 


} 


/* */ 

get_dvect (xp , yp , ncs , nseg , n , if 1 , brkpt , dvect_in , dvect_out ) 

float xp [ ] [MAXCS ] , yp [ ] [MAXCS ] , dvect_in [ ] [MAXCS ] [ 2 ] , dvect_out [ ] [MAXCS ] [ 2 ] ; 
int ncs, nseg, n[] [MAXCS] ,ifl[] [MAXCS] , brkpt [] [MAXCS] ; 


float a[4] ,b[4] , ct,st,v[2] ,len,xel,xe2,yel,ye2; 
int cs, I,seg,nl,n2, ibase; 


{ 
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for (ca - 0 ; cs < ncs ; cs++) 


/* loop over all cross sections */ 


/* loop over all segment end points 
for (seg - 0 ; seg <■ nseg ; seg++) { 

nl « n[seg][cs]; /* # of points in prior segment 

n2 » n[seg + 1] [cs] ; /* # of points in following segment 

ibase - brkpt[seg] [cs] ; /* subscript of the breakpoint 

/* CASE: continuous derivative */ 

/* neither curve linear */ 

if ((nl > 2) && (n2 > 2) && (ifl[ibase] [cs] — l) ) { 

a[0] - xp[ibase - 1] [cs] - xp[ibase - 2][cs]; 
a[l] - xp[ibase] [cs] - xp[ibase - l][cs]; 
a[2] * xp[ibase + 1] [cs] - xp[ibase] [cs] ; 
a[3] = xp[ibase + 2] [cs] - xp[ibase + 1] [cs] ; 

b[0] ■ yp[ibase - 1] [cs] - yp[ibase - 2][cs]; 

b [l] “ ypfibase] [cs] - yp[ibase - 1] [cs] ; 
b[2 ] ■ yp[ibase + 1] [cs] - yp[ibase] [cs ] ; 

b[3] ■ yp[ibase + 2][cs] - yp[ibase + 1] [cs] ; 

tancalc(a,b / &ct / &st) ; 

dvect_in[seg] [cs] [0] * -ct; 
dvect_in[seg] [cs] [1] ■ -st; 
dvect_out [ seg ] [ cs ] [ 0 ] * ct ; 
dvect_out [ seg ] [ cs ] [ 1 ] = st ; 

> 

- /* CASE: continuous derivative */ 

/* prior curve linear */ 

else if ((nl — 2) && ( ifl[ ibase ] [cs] *=1)) { 

v[0] - xp[ibase - 1] [cs] - xp[ibase] [cs] ; 

v [ X ] » ypfibase - 1] [cs] - yp [ ibase] [cs] ; 

len » ln_vect(v) ; 

ct = v [ 0 ] / len; 

st * v[l] / len; 

dvect_in[seg] [cs] [0] *=ct; 

dvect_in[seg] [cs] [1] » st; 

dvect_out [ seg ][cs][0] = -ct; 

dvect_out [seg] [cs] [1] ■ -st; 


/* CASE: continuous derivative */ 
/* following curve linear */ 
else if ( (n2 — 2) && (ifl[ibase] [cs] ===== 1) ) { 

v[0] = xp[ibase + 1] [cs] - xp[ ibase] [cs] ; 
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v[l] ■ ypfibase + 1] [cs] - yp[ibase] [cs] ; 
len - ln_vect(v) ; 
ct ** v[0] / len; 
st ■ v[l] / len; 
dvect_in[seg] [cs] [0] * -ct; 
dvect_in[seg] [cs] [1] * -st; 
dveqt_out [ seg ] [ cs ] [ 0 ] * ct; 
dvect_out[seg] [cs] [1] » st; 


/* all remaining CASES have */ 

/* discontinuous derivative */ 

else { 

/* incoming derivative calculation */ 
/* prior curve linear */ 

if (nl — 2) { 

v[0] = xp[ibase - l][cs] - xp[ ibase] [cs] ; 
v[l] - yp[ibase - l][cs] - yp[ibase] [cs] ; 
len = ln_vect(v) ; 

dvect_in[seg] [cs] [0] * v[0] / len; 
dvect_in[seg] [cs] [1] - v[l] / len; 


/* incomming derivative calculation */ 
/* prior curve non-linear */ 

else if (seg != 0) { 

xel ■ extrap ( xp [ ibase ][cs],xp[ ibase - 1] [cs] /Xp[ibase - 2][cs]) 

xe2 » extrap (xel ,xp[ ibase] [cs] ,xp[ ibase - l][cs]); 

yel = extrap (yp[ ibase ] [cs] ,yp[ ibase - 1] [cs] ,yp[ibase - 2][cs]) 

ye2 - extrap ( yel ,yp[ ibase ] [cs] ,yp[ ibase - l][cs]); 

a[0] » xp[ibase - l][cs] - xp[ibase - 2][cs]; 

a[l] - xp [ ibase ] [cs] - xp[ibase - 1] [cs] ; 

a[2] ■ xel - xp[ ibase] [cs] ; 

a[3] * xe2 - xel; 

• k[0] " yp[ibase - 1] [cs] - yp[ibase - 2][cs]; 
b[l] ■ yp [ ibase ] [cs] - yp[ibase - l][cs]; 
b [2] * yel - yp[ibase] [cs] ; 
to C 3 3 « ye2 - yel; 

tancalc(a,b / &ct / &st) ; 

dvect_in[seg] [cs] [0] = -ct; 
dvect_in[seg] [cs] [l] ■ -st; 

• > 


if (n2 == 2) { 


/* outgoing derivative calculation */ 
/* following curve linear */ 
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v[0] » xp[ibase + l][cs] - xp[ ibase] [cs] ; 
v[l] “ yp[ibase + l][cs] - yp[ibase] [cs] ; 
len - ln_vect(v) ; 

dvect_out[seg] [cs] [0] = v[0] / len? 
dvect_out[seg] [cs] [1] = v[l] / len; 


else if (seg 1= nseg) 


/* outgoing derivative calculation 
/* following curve non-linear 


xel - extrap (xp[ ibase] [cs] ,xp[ibase+l] [cs] ,xp[ibase+2] [cs] ) ; 

xe2 - extrap (xel, xp[ ibase] [cs] ,xp[ibase+l] [cs] ) ; 

y,el - extrap (yp[ ibase] [cs] ,yp[ibase+l] [cs] ,yp[ibase+2] [cs] ) ; 

y®2 - extrap (yel,yp[ ibase] [cs] ,yp[ibase+l] [cs] ) ; 

a[0] ■ xel - xe2? 

a[l] - xp[ibase] [cs] - xel; 

a[2] - xp[ ibase + l][cs] - xp[ibase] [cs] ? 

a[3] - xp[ ibase + 2][cs] - xp[ibase + l] [cs] ; 

b[0] ■ yel - ye2; 

b[l] * yp[ibase] [cs] - yel; 

b[2] - yp[ibase + l][cs] - yp[ibase] [cs] ; 

b[3] = yp[ibase + 2][cs] - yp[ibase + 1] [cs] ; 

tancalc ( a , b , &ct , &st ) ; 

dvect_out [seg] [cs] [0] ■ ct? 
dvect__out[seg] [cs] [1] ■ st; 


} 

) /* end loop over break points */ 


ct = dvect_out[0] [cs] [0] ; 
st * dvect_out[0] [cs] [1] ; 
setscl(&ct, &st) ; 
dvect_out [ 0 ] [ cs ] [ 0 ] ■ ct ; 
dvect~out[0] [cs] [1] = st; 

ct = dvect_in[nseg] [cs] [0] ; 
st = dvect_in[nseg] [cs] [l] ? 
setscl(&ct, &st) ; 

♦ dvect_in[nseg] [cs] [0] = ct; 
dvect_in[nseg] [cs] [1] = st; 

} /* end loop over cross sections * 


/* reset the first and last */ 

/* derivative values as required */ 
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/* */ 

bz_minerr 1 (x , y , do , dl , ncs , n , segno , lambda_old , lambda) 

float x[] [MAXCS], y[] [MAXCS] ,d0[] [2],dl[] [2] ,lambda_old[] , lambda [) [2] ; 
int ncs , n [ ] [MAXCS ] , segno ; 

{ float dist, s [MAXCS] [5] ,scomp [MAXCS] [5] ,sx2m3 [MAXCS] ,sx2pl [MAXCS] ,11,12,13 
rl , r2 , r3 , r4 , r5 , r6 , r , totl , tot2 , totald , 
a[MAXCS] , b [MAXCS] ,C[MAXCS],d [MAXCS] , f [MAXCS] ; 
int i,cs,p,n_cs,nml,nm2 ; 


/* loop over all cross sections */ 


/* number of points in the cross */ 
/* section */ 

/* straight line => no error */ 


/* error possible in all other cases 
/* find chord lengths of data */ 
/* positions & related functions */ 
/* required for matrix set up */ 
for (i - 0 ; i < nml ; i++) ( 

- dist * sqrt(sqr(x[i+l] [cs] - x[i][cs]) 

+ sqr(y[i+l] [cs] - y[i][cs])); 
totald +* dist; 
s[i] [0] ■ totald; 

) 

for (i ■ 0 ; i < nm2 ; i++) { 
s[i] [0] /= totald; 
scomp[i][0] * 1 - s[i] [0] ; 
sx2m3 [i] = 3 - 2 * s[i] [0] ; 

sx2pl[i] - 1 + 2 * s[i] [0] ; 

for (p ■ 1 ; p < 4 ; p++) { 

s[i][p] - s[i][p - 1] * s[i] [0] ; 
scomp[i][p] = scomp[i] [p-1] * scomp[i][0]; 


for (cs - 0 ;. cs < ncs ; cs++) { 
totald « 0; 

n_cs ■ n[segno + l][cs]; 
nml ■ n_cs -1 ; 
nm2 - nml - 1; 


switch (n_cs) { 


case 2: 
break; 

default: 
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} 


11 - o; 

12 * 0 ; 
13 - 0; 
rl - o; 
r2 * 0; 
r3 » 0; 
r4 - 0; 
r5 ** 0; 
r6 ■ 0; 


/* set up the least squares matrix */ 
/* with information from this */ 
/* cross section */ 


for (i * 0 ; i < nm2 ; i++) { 

11 +- scomp[i] [3 ] * s[i] [1] ; 

12 +- scomp[i] [2] * s[i] [2] ; 

13 +» scomp[i] [1] * s[i] [3] ? 

rl +- scomp[i] [1] * s[i] [0] * (x[i + 1] [cs] * d0[cs][0] 

+ y[i + 1] [cs] * d0[cs][l]); 
r2 +■ scomp[i][3] * s[i][0] * sx2pl[i]; 

r3 +- scomp [ i ] [ 1 ] * s [ i ] [ 2 ] * sx2m3 [ i ] ; 

r4 +* scomp[i][0] * s[i][l] * (x[i + 1] [cs] * dl[cs][0] 

+ y[i + l][cs] * dl[cs] [1] ) ; 
r5 +«■ scomp[i][2] * s[i][l] * sx2pl[i]; 

r6 +- scomp [i][0] * s[i][3] * sx2m3[i]; 


a[cs] ■ 3 * 12 * (d0[cs] [0] * dl[cs][0] + d0[cs][l] * dl[cs][l]) 

* lambda_old[cs] ; 

b[cs] - 3 * 13 * (sqr (dl[cs] [0] ) + sqr (dl[cs] [1] ) ) ; 
c [cs] ■ r4 - r5 * (x[0] [cs] * dl[cs][0] + y[0] [cs] * dl[cs][l]) 

- r6 * (x[nml] [cs] * dl[cs] [0] + y[nml] [cs] * dl[csj[l]); 
d[cs] * 3 * 11 * (sqr(d0[cs] [0] ) + sqr (d0[cs] [1] ) ) 

* sqr (lambda_old[cs] ) ; 

f[cs] - (rl - r2 * (x[0] [cs] * d0[cs][0] + y[0][cs] * d0[cs][l]) 

- r3 * (x[nml][cs] * d0[cs][0] + y[nml][ cs i * d0[cs][l])) 

* lambda_old[cs] ; 

break; 

} 


} /* end loop over cross sections */ 

m 

/* solve the first curve segment */ 
if (segno —*■ 0) /* no prior "r" constraint */ 

«solvel(a,b,c,d, f ,ncs, n,x,y, lambda) ; 

/* solve other curve segments * 

else /* use common "r" at each start point * 

solve2 ( a , b , c , d , f , ncs , segno , n , lambda_old , x , y , lambda ) ; 
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/* 

float ln_vect(v) 
float v[2] ; 

{ float len; 

len = sqrt(sqr(v[0] ) + sqr(v[l])); 
return (len) ; 


/* 

setscl(ct,st) 
float *ct-, *st; 
{ 


if (fabs(*ct) < .05) { 

*ct - 0; 

*st - 1; 

) 

else if (fabs(*st) < .05) ( 
*ct - 1; 

*st = 0; 

} 


> 


/* 

float extrap ( a, b f c) 
float a,b,c; 

{ float d; 

d * 3 * (a-b) + c; 


return d; 
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} 


tancalc (a,b,dcos,dsin) 

float a[4] ,b[4] ,*dcos f *dsin; 

{ float a0,b0,tl,t2,w2,w3? 

w2 - fabs(a[2] * b[3] - a[3] * b[2])y 

w3 - fabs (a[0] * b[l] - a[l] * b[0]); 

aO « w2 * a[l] + w3 * a[2]; 

bO - w2 * b[l] + w3 * b[2] ; 

tl - sqrt(a0 * aO + bO * bO) ; 


/* If the curve is a straight line, then */ 
/* b[i]/a[i] = b[i+l]/a[i+l] = tan(a) V 
/* for all i = 0,1,2 */ 
/* and thus aO = bo * 0 */ 
/* treat this case separately. */ 


if ((aO — 0) && (bO — 0)) { 

t2 - sqrt(sqr(a[0] ) + sqr(b[0])); 
*dcos * a[0] / t2; 

*dsin - b[0] / t2 ; 

} 

else { 

*dcos * aO / tl? 

*dsin ■ bO / tl; 


> 


/* */ 

sep_seg (xp , yp , dvect_in , dvect_out , brkpt , segno , ncs , x , y , do , dl ) 

float xp[] [MAXCS] ,yp[] [MAXCS] , dvect_in [ ] [MAXCS ] [2] ,dvect_out[ ] [MAXCS] [2] , 
X[] [MAXCS] ,y[] [MAXCS] ,d0[] [2] ,dl[] [2] ; 
int brkpt [ ] [MAXCS ] , segno , ncs ? 

{ int cs, istart, iend, i, j ? 

/* loop over all cross sections */ 
for (cs - 0 ? cs < ncs ? cs++) { 
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/* subscripts for the start and end */ 
/* of the current curve segment in */ 
/* this cross section */ 

istart - brkpt [ segno ] [cs] ; 
iend * brkpt [segno + 1] [cs] ; 
j - 0? 


/* loop over points of this segment */ 
/* copy to "x" & "y M */ 

for (i * istart ; i <= iend ; i++) { 

x[j] [cs] = xp[i] [cs] ; 
y[j][cs] » yp[i] [cs] ; 
j++; 

} /* end loop over points */ 

/* copy derivatives for this segment */ 
d0[cs][0] * dvect_out [segno] [cs] [0] ? 
d0[cs] [1] = dvect_out[ segno] [cs] [1] ; 
dl[cs][0] = dvect_in[ segno + l][cs][0]; 
dl[cs] [1] - dvect_in[ segno + l][cs][l]; 

} /* end loop ove cross sections */ 


} 


/* */ 

solve 1 (a/bfC/dff, ncs ,n,x,y, lambda) 

float a[],b[],c[] / d[] / f[] / x[] [MAXCS],y[] [MAXCS] , lambda [] [2] ; 
int ncs,n[] [MAXCS] ? 

{ int cs ; 

float det,v[2]; 


/* loop over cross sections */ 
for (cs m o ; cs < ncs ; cs++) { 




/* CASE: non-linear curve segment */ 
/* solve the matrix */ 

if (n[l] [cs] != 2) { 

det = d[cs] * b[cs] - sqr(a[cs]); 

lambda [cs] [0] * (f[cs] * b[cs] - c[cs] * a[cs]) / det; 
lambda[cs] [1] « (c[cs] * d[cs] - f[cs] * a[cs]) / det; 


else { 

v[0] - x[l][cs] 
v[l] = y[l] [cs] 


/* CASE: linear segment */ 

x[0] [cs] ; 
y[0] [cs] ; 
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lambda [cs] [0] - ln_vect(v) / 3.; 
lambda [cs] [1] - lambda[cs] [0 ] ; 

} 

) /* end loop over cross sections */ 

} 


/* */ 

solve2 ( a , b , c , d , f , ncs , segno , n , lambda_old , x , y , lambda ) 

float a[] ,b[] ,c[] ,d[],f [] , lambda_old[ ] ,x[] [MAXCS], y[] [MAXCS] ,lambda[] [2] ; 
int ncs , segno , n [ ] [MAXCS ] ; 

{ float totl,tot2,r,v[2] ; 
int cs,segpl; 

totl - 0; 
tot2 * 0; 

segpl - segno + 1; 

/* solve first for the required common */ 

/* ratio */ 

/* loop over all cross sections */ 
for (cs - 0 ; cs < ncs ; cs++) { 

if (n[segpl] [cs] 1= 2) { 

totl +- f[cs] - (a[cs] *c[cs]) / b[cs]; 
tot2 += d[csj - sqr(a[cs]) / b[cs] ; 

) 

} /* end loop over cross sections */ 

r » totl / tot2y 
* 

/* use "r M to generate the first control */ 
/* point, then solve for the other */ 

/* loop over all cross sections */ 

for (cs - 0 ; cs < ncs ? cs++) { 

/* first control point */ 
lambda [cs] [0] ■ r * lambda_old[cs] ; 

' /* second control point */ 

/* CASE: non-linear segment */ 

if (n[segpl] [cs] 1* 2) 

lambda[cs] [1] ■ (c[cs] - (a[cs] * r) ) / b[cs]; 

/* CASE: linear segment */ 
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v[0] - x[l][cs] - x[0] [cs] ; 
v[l] - y[l][cs] - y[0][cs]; 
lambda [cs] [1] * In vect(v) / 3.; 

> 

} /* end loop over cross sections */ 


} 


/* V 

get_cpts ( xp , yp / zp , do , dl , brkpt , segno , ncs , lambda , cpts ) 

float xp[] [MAXCS] ,yp[] [MAXCS] , zp[] ,d0[] [2] ,dl[] [2] , lambda [] [2] , 
cpts [ ] [MAXCS ] [4] [3] ; 
int brkpt [ ] [MAXCS ] , segno , ncs ; 

{ int cs , istart , iend , i ; 


/* store the control points of this */ 

/* segment for each cross section */ 

/* loop over each cross section */ 
for (cs ■ 0 ? cs < ncs ; cs++) { 

istart « brkpt [ segno ] [cs] ; 
iend = brkpt [segno + 1] [cs] ? 

cpts [segno] [cs] [0] [0] * xp [ istart ] [cs] ; 
cpts[segno] [cs] [0] [1] ■ yp [ istart ] [cs] ; 

cpts[segno] [cs] [1] [0] = xp [ istart ] [cs] + lambda[cs] [0] * d0[cs][0] 
cpts[segno] [cs] [1] [1] * yp [ istart ] [cs] + lambdafcs] [0] * d0[cs][l] 
cpts [ segno] [cs] [2] [0] » xp[iend][cs] + lambda[cs] [1] * dl[cs][0]; 
cpts[segno] [cs] [2] [1] * yp[iend][cs] + lambda[cs] [1] * dl[cs][l]; 
cpts[segno] [cs] [3] [0] * xp[iend] [cs] ; 
cpts[segno] [cs] [3] [1] * yp[iend] [cs] ; 

for (i = 0 t i < 4 ; i++) 

cpts[segno] [cs] [i] [2] - zp[cs]; 

) 

/* end loop over cross sections */ 


% 


} 


